! *****************************************************************************
! Setup Grid and Paritcle Parameters
! *****************************************************************************

module Setup

    use mpi
    use hdf5

    implicit none

! ******************************************************
! Physical Constants in cgs Unit
! ******************************************************

    real(kind=8), parameter :: pi = 3.14159d0  ! Circular Constant
    real(kind=8), parameter :: kB = 1.3807d-16 ! Boltzmann Constant [erg/K]
    real(kind=8), parameter :: c = 2.9979d10   ! Light Speed [cm/s]
    real(kind=8), parameter :: e = 4.8032d-10  ! Elementary Charge [statcoul]
    real(kind=8), parameter :: me = 9.1094d-28 ! Electron Mass[g]
    real(kind=8), parameter :: mp = 1.6726d-24 ! Proton Mass [g]

! ******************************************************
! Internal Parameters
! ******************************************************

! Define Particle Linked List
    type PseudoParticle
        integer(kind=4) :: Z, A
        real(kind=8) :: x, y
        real(kind=8) :: vx, vy, vz
        real(kind=8) :: weight
        type(PseudoParticle), pointer :: next
    end type PseudoParticle

    type(PseudoParticle), pointer :: IonListHead=>null()
    type(PseudoParticle), pointer :: InjectorListHead=>null()
    integer(kind=4) :: LocalIons  ! Particle Number per Processor

! Define Fluid Particle
    type FluidParticle
        integer(kind=4) :: Z, A
        integer(kind=4) :: ppc
        real(kind=8) :: ux, uy, uz
        real(kind=8) :: Ti, Te
        real(kind=8) :: ni
    end type FluidParticle

    real(kind=8) :: A_avg, Z_avg

! Define Grid Variables
    real(kind=8), allocatable :: x(:,:) ! Grid Coordinate in X
    real(kind=8), allocatable :: y(:,:) ! Grid Coordinate in Y
    real(kind=8), allocatable :: ne(:,:) ! Electron Number Density
    real(kind=8), allocatable :: ni(:,:) ! ion Number Density
    real(kind=8), allocatable :: ne_prev(:,:) ! Electron Number Density at last time step
    real(kind=8), allocatable :: ne_temp(:,:) ! Electron Number Density at time step N+1/2
    real(kind=8), allocatable :: Te(:,:) ! Electron Temperature
    real(kind=8), allocatable :: ux(:,:) ! Fluid Velocity in X
    real(kind=8), allocatable :: uy(:,:) ! Fluid Velocity in Y
    real(kind=8), allocatable :: uz(:,:) ! Fluid Velocity in Z
    real(kind=8), allocatable :: pe(:,:) ! Electron Pressure
    real(kind=8), allocatable :: Ex(:,:) ! Electric Field in X
    real(kind=8), allocatable :: Ey(:,:) ! Electric Field in Y
    real(kind=8), allocatable :: Ez(:,:) ! Electric Field in Z
    real(kind=8), allocatable :: Exg(:,:) ! Electric Field in X
    real(kind=8), allocatable :: Eyg(:,:) ! Electric Field in Y
    real(kind=8), allocatable :: Bx0(:,:) ! Electric Field in X to calculate Bf
    real(kind=8), allocatable :: By0(:,:) ! Electric Field in Y to calculate Bf
    real(kind=8), allocatable :: Bz0(:,:) ! Electric Field in Z to calculate Bf
    real(kind=8), allocatable :: Bx(:,:) ! Total Magnetic Field in X on half grid
    real(kind=8), allocatable :: By(:,:) ! Total Magnetic Field in Y on half grid
    real(kind=8), allocatable :: Bz(:,:) ! Total Magnetic Field in Z on half grid
    real(kind=8), allocatable :: Bxg(:,:) ! Total Magnetic Field in X
    real(kind=8), allocatable :: Byg(:,:) ! Total Magnetic Field in Y
    real(kind=8), allocatable :: Bzg(:,:) ! Total Magnetic Field in Z
    real(kind=8), allocatable :: Bwx(:,:) ! Ion-Weibel Magnetic Field in X
    real(kind=8), allocatable :: Bwy(:,:) ! Ion-Weibel Magnetic Field in Y
    real(kind=8), allocatable :: Bwz(:,:) ! Ion-Weibel Magnetic Field in Z
    real(kind=8), allocatable :: Bfx(:,:) ! Faraday Magnetic Field in X
    real(kind=8), allocatable :: Bfy(:,:) ! Faraday Magnetic Field in Y
    real(kind=8), allocatable :: Bfz(:,:) ! Faraday Magnetic Field in Z

    logical :: UseConv, UseCAMCL, UseCNPPC

! Define EOS Variables
    ! real(kind=8) :: Gamma = 1.666667d0 ! Ratio of Specific Heat
    real(kind=8) :: Gamma = 5.0/3.0 ! Ratio of Specific Heat
    real(kind=8), allocatable :: ne0(:,:) ! Electron Number Density
    real(kind=8), allocatable :: Te0(:,:) ! Electron Temperature
    real(kind=8), allocatable :: pe0(:,:) ! Electron Pressure

! ******************************************************
! Parameters in Control Block
! ******************************************************

! Running Time
    real(kind=8) :: TotalTime
    real(kind=8) :: dt
    real(kind=8) :: RunTime
    real(kind=8) :: StartTime
    real(kind=8) :: EndTime
    real(kind=8) :: Nt_unit
    integer(kind=4) :: nt_b     ! sub-cycle for B field

! Size of Simulation Box
    real(kind=8) :: xMinBox, xMaxBox, yMinBox, yMaxBox, CxDebye, CyDebye
    real(kind=8) :: dx, dy
    real(kind=8) :: Nx_unit, Ny_unit
    integer(kind=4) :: NxBox, NyBox

! MPI Parameters
    integer(kind=4) :: ppx, ppy
    integer(kind=4) :: TotalProcessors ! Total Number of Processors
    integer(kind=4) :: myid ! Process Number
    integer(kind=4) :: status(mpi_status_size), ierr
    integer(kind=4) :: ppi, ppj
    integer(kind=4) :: Nx, Ny, ncx, ncy, inx, iny
    real(kind=8) :: xMin, xMax, yMin, yMax
    type(PseudoParticle), pointer :: MPIListHeadLeft=>null()
    type(PseudoParticle), pointer :: MPIListHeadRight=>null()
    type(PseudoParticle), pointer :: MPIListHeadUp=>null()
    type(PseudoParticle), pointer :: MPIListHeadDown=>null()
    type(PseudoParticle), pointer :: MPIInjectorListHeadLeft=>null()
    type(PseudoParticle), pointer :: MPIInjectorListHeadRight=>null()
    type(PseudoParticle), pointer :: MPIInjectorListHeadUp=>null()
    type(PseudoParticle), pointer :: MPIInjectorListHeadDown=>null()
    integer(kind=4) :: MPIIonsLeft  ! Size of MPIListHeadLeft
    integer(kind=4) :: MPIIonsRight ! Size of MPIListHeadRight
    integer(kind=4) :: MPIIonsUp    ! Size of MPIListHeadUp
    integer(kind=4) :: MPIIonsDown  ! Size of MPIListHeadDown

! Dimensionless Parameters  use [cgs]
    real(kind=8) :: L_unit, time_unit, v_unit 
    real(kind=8) :: E_unit, B_unit, M_unit
    real(kind=8) :: n_unit ! = M_unit*c**2/4/pi/e**2/L_unit**2
    real(kind=8) :: p_unit ! = M_unit*c**2/4/pi/e**2/time_unit**2
    real(kind=8) :: T_unit ! = M_unit*v_unit**2/kB/11605 ![eV]
    real(kind=8) :: TeCharacter ! [eV]
    real(kind=8) :: neCharacter ! [eV]
    real(kind=8) :: Omega_p ! = sqrt(4*pi*(Z*e)**2*neCharacter/M_unit)
    real(kind=8) :: Omega_g ! = B_unit*Z*e/M_unit/c
    
    real(kind=8) :: lambda_De ! = sqrt(kB*TeCharacter/4/pi/neCharacter/e**2)
    real(kind=8) :: SkinDepth ! = c/Omega_p (inertial length)
    real(kind=8) :: v_th ! = sqrt(kB*TeCharacter/M_unit)
    real(kind=8) :: v_Alfen ! = B_unit/sqrt(4*pi*neCharacter*M_unit)
    real(kind=8) :: Cs ! = sqrt(Gamma*Z*kB*TeCharacter/M_unit)
    real(kind=8) :: gyroradius ! = v_th/Omega_g
    real(kind=8) :: electronTrap

! ******************************************************
! Parameters in Boundary Block
! ******************************************************
    logical :: UseBoundary
    integer(kind=4) :: xMinBoundary
    integer(kind=4) :: xMaxBoundary
    integer(kind=4) :: yMinBoundary
    integer(kind=4) :: yMaxBoundary
! ******************************************************
! Parameters in Plasma Block
! ******************************************************
    logical :: UsePlasma
    integer(kind=4) :: NSpecies

! ******************************************************
! Parameters in Field Block
! ******************************************************
    logical :: UseField

! ******************************************************
! Parameters in Injector Block
! ******************************************************
    logical :: UseInjector
    integer(kind=4) :: NInjector
    logical :: xMinInjector, xMaxInjector
    logical :: yMinInjector, yMaxInjector
    real(kind=8) :: xInjectMin,xInjectMax,yInjectMin,yInjectMax

! ******************************************************
! Parameters in Physics Block
! ******************************************************
    logical :: UseLorentz
    logical :: UseWeibel
    logical :: UseFaraday
    logical :: UseOhm
    logical :: UseICS

! ******************************************************
! Parameters in Output Block
! ******************************************************
    integer(kind=4) :: DumpNumber
    real(kind=8) :: DumpPeriod
    logical :: neDump       ! Density
    logical :: niDump       ! Density
    logical :: uxDump       ! Velocity in X
    logical :: uyDump       ! Velocity in Y
    logical :: uzDump
    logical :: peDump       ! Electron Pressure
    logical :: ExDump       ! Electric Field in X
    logical :: EyDump       ! Electric Field in Y
    logical :: EzDump       ! Electric Field in Z
    logical :: BxDump       ! Magnetic Field in X
    logical :: ByDump       ! Magnetic Field in Y
    logical :: BzDump       ! Magnetic Field in Z
    logical :: BwxDump      ! Ion-Weibel Magnetic Field in X
    logical :: BwyDump      ! Ion-Weibel Magnetic Field in Y
    logical :: BwzDump      ! Ion-Weibel Magnetic Field in Z
    logical :: IonsDump     ! Particle Ions

    integer, parameter :: ndim = 2
    integer(HSIZE_T), dimension(ndim) :: dimsf
    integer(HID_T) :: file_id       ! File identifier
    integer(HID_T) :: dset_id       ! Dataset identifier
    integer(HID_T) :: dataspace     ! Dataspace identifier
    integer(HID_T) :: memspace      ! memspace identifier
    integer(HID_T) :: plist_id      ! Property list identifier 
    integer(HSIZE_T), dimension(ndim)  :: datasize
    integer(HSIZE_T), dimension(ndim)  :: offset
 
    integer :: mpierror,error       ! MPI error flag
    integer :: mpi_size, mpi_rank

! ******************************************************
! Parameters in Dimensionless Block
! ******************************************************
    logical :: UseDimensionless

! ******************************************************
! Parameters in Time Block
! ******************************************************
    logical :: UseTime

contains

!!!!!! Functions for Ion Linked-List !!!!!!!!!!!!!!!!!!!!!!
    subroutine createHeadOf(p)
        implicit none
        type(PseudoParticle), pointer :: p, q

        allocate(q)
        q%Z = 0
        q%A = 0
        q%x = 0.0d0
        q%y = 0.0d0
        q%vx = 0.0d0
        q%vy = 0.0d0
        q%vz = 0.0d0
        q%weight = 0.0d0
        nullify(q%next)
        p => q
    end subroutine createHeadOf

    subroutine isInside(Ion,IonIsInside)
        implicit none
        type(PseudoParticle), pointer :: Ion
        logical :: IonIsInside

        if(Ion%x<xMin)then
            call insertIon(MPIListHeadLeft,Ion)
            IonIsInside = .false.
        elseif(Ion%x>xMax)then
            call insertIon(MPIListHeadRight,Ion)
            IonIsInside = .false.
        elseif(Ion%y>yMax)then
            call insertIon(MPIListHeadUp,Ion)
            IonIsInside = .false.
        elseif(Ion%y<yMin)then
            call insertIon(MPIListHeadDown,Ion)
            IonIsInside = .false.
        else
            IonIsInside = .true.
        end if
    end subroutine isInside

    subroutine isInInjector(Ion,InjectorIsInside)
        implicit none
        type(PseudoParticle), pointer :: Ion
        logical :: InjectorIsInside

        if(Ion%x<xInjectMin)then
            call insertIon(MPIInjectorListHeadLeft,Ion)
            InjectorIsInside = .false.
        elseif(Ion%x>xInjectMax)then
            call insertIon(MPIInjectorListHeadRight,Ion)
            InjectorIsInside = .false.
        elseif(Ion%y>yInjectMax)then
            call insertIon(MPIInjectorListHeadUp,Ion)
            InjectorIsInside = .false.
        elseif(Ion%y<yInjectMin)then
            call insertIon(MPIInjectorListHeadDown,Ion)
            InjectorIsInside = .false.
        else
            InjectorIsInside = .true.
        end if
    end subroutine isInInjector

    subroutine isInInjector2(Ion,InjectorIsInside)
        implicit none
        type(PseudoParticle), pointer :: Ion
        logical :: InjectorIsInside

        if(Ion%x<=xInjectMax .and. Ion%x>=xInjectMin &
        .and. Ion%y<=yInjectMax .and. Ion%y>=yInjectMin)then
            InjectorIsInside = .true.
        else
            InjectorIsInside = .false.
        end if
    end subroutine isInInjector2

    subroutine deleteIon(ListHead,PrevIon)
        implicit none
        type(PseudoParticle), pointer :: ListHead, PrevIon
        type(PseudoParticle), pointer :: RemovedIon

        RemovedIon => PrevIon%next
        PrevIon%next => RemovedIon%next
        deallocate(RemovedIon)
        if(associated(ListHead,IonListHead)) LocalIons = LocalIons-1
        if(associated(ListHead,MPIListHeadLeft)) MPIIonsLeft = MPIIonsLeft-1
        if(associated(ListHead,MPIListHeadRight)) MPIIonsRight = MPIIonsRight-1
        if(associated(ListHead,MPIListHeadUp)) MPIIonsUp = MPIIonsUp-1
        if(associated(ListHead,MPIListHeadDown)) MPIIonsDown = MPIIonsDown-1
    end subroutine deleteIon

    subroutine deleteList(ListHead)
        implicit none
        type(PseudoParticle), pointer :: ListHead
        type(PseudoParticle), pointer :: p, q

        q => ListHead
        p => ListHead%next
        do while(associated(p))
            call deleteIon(ListHead,q)
            p => q%next
        end do
    end subroutine deleteList

    subroutine insertIon(ListHead,Ion)
        implicit none
        type(PseudoParticle), pointer :: ListHead, Ion
        type(PseudoParticle), pointer :: p

        allocate(p)
        p%A = Ion%A
        p%Z = Ion%Z
        p%x = Ion%x
        p%y = Ion%y
        p%vx = Ion%vx
        p%vy = Ion%vy
        p%vz = Ion%vz
        p%weight = Ion%weight
        
        p%next => ListHead%next
        ListHead%next => p
        if(associated(ListHead,IonListHead)) LocalIons = LocalIons+1
        if(associated(ListHead,MPIListHeadLeft)) MPIIonsLeft = MPIIonsLeft+1
        if(associated(ListHead,MPIListHeadRight)) MPIIonsRight = MPIIonsRight+1
        if(associated(ListHead,MPIListHeadUp)) MPIIonsUp = MPIIonsUp+1
        if(associated(ListHead,MPIListHeadDown)) MPIIonsDown = MPIIonsDown+1
    end subroutine insertIon

    subroutine copyIon(ListHead,Ion)
        implicit none
        type(PseudoParticle), pointer :: ListHead, Ion
        type(PseudoParticle), pointer :: p

        allocate(p)
        p%A = Ion%A
        p%Z = Ion%Z
        p%x = Ion%x
        p%y = Ion%y
        p%vx = Ion%vx
        p%vy = Ion%vy
        p%vz = Ion%vz
        p%weight = Ion%weight
        
        p%next => ListHead%next
        ListHead%next => p

    end subroutine copyIon

    subroutine sendIons(obj,ListHead)
        implicit none
        integer(kind=4) :: obj
        type(PseudoParticle), pointer :: ListHead
        integer(kind=4) :: num, i
        real(kind=8), allocatable :: Temp(:,:)
        type(PseudoParticle), pointer :: p, q

        if(associated(ListHead,MPIListHeadLeft)) num = MPIIonsLeft
        if(associated(ListHead,MPIListHeadRight)) num = MPIIonsRight
        if(associated(ListHead,MPIListHeadUp)) num = MPIIonsUp
        if(associated(ListHead,MPIListHeadDown)) num = MPIIonsDown

        call MPI_SEND(num,1,mpi_integer,obj,1,mpi_comm_world,ierr)
        if(num>0)then
            allocate(Temp(num,8))
            q => ListHead
            p => ListHead%next
            do i=1,num
                Temp(i,1) = real(p%A,8)
                Temp(i,2) = real(p%Z,8)
                Temp(i,3) = p%x
                Temp(i,4) = p%y
                Temp(i,5) = p%vx
                Temp(i,6) = p%vy
                Temp(i,7) = p%vz
                Temp(i,8) = p%weight
                call deleteIon(ListHead,q)
                p => q%next
            end do
            call MPI_SEND(Temp,num*8,mpi_real8,obj,1,mpi_comm_world,ierr)
            deallocate(Temp)
        end if
    end subroutine sendIons

    subroutine recvIons(obj,ListHead)
        implicit none
        integer(kind=4) :: obj
        type(PseudoParticle), pointer :: ListHead
        integer(kind=4) :: num, i
        real(kind=8), allocatable :: Temp(:,:)
        type(PseudoParticle), pointer :: p
        logical :: IonIsInside

        call MPI_RECV(num,1,mpi_integer,obj,1,mpi_comm_world,status,ierr)
        if(num>0)then
            allocate(Temp(num,8))
            call MPI_RECV(Temp,num*8,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            allocate(p)
            do i=1,num
                p%A = int(Temp(i,1),4)
                p%Z = int(Temp(i,2),4)
                p%x = Temp(i,3)
                p%y = Temp(i,4)
                p%vx = Temp(i,5)
                p%vy = Temp(i,6)
                p%vz = Temp(i,7)
                p%weight = Temp(i,8)
                call isInside(p,IonIsInside)
                if(IonIsInside) call insertIon(IonListHead,p)
            end do
            deallocate(Temp)
            deallocate(p)
        end if
    end subroutine recvIons

    subroutine moveListA2B(ListHeadA,ListHeadB)
        implicit none
        type(PseudoParticle), pointer :: ListHeadA, ListHeadB
        type(PseudoParticle), pointer :: p, q

        q => ListHeadA
        p => ListHeadA%next
        do while(associated(p))
            call insertIon(ListHeadB,p)
            call deleteIon(ListHeadA,q)
            p => q%next
        end do
    end subroutine moveListA2B

    subroutine copyListA2B(ListHeadA,ListHeadB)
        implicit none
        type(PseudoParticle), pointer :: ListHeadA, ListHeadB
        type(PseudoParticle), pointer :: p, q

        q => ListHeadA
        p => ListHeadA%next
        do while(associated(p))
            call insertIon(ListHeadB,p)
            q => p
            p => p%next
        end do
    end subroutine copyListA2B

!!!!!!! The Mapping Functions for Particle and Grid !!!!!!!
    function mapParticle2Grid(rp,U,weight,Z)
        implicit none
        real(kind=8) :: rp(2), weight, U
        integer(kind=4) :: Z
        real(kind=8) :: mapParticle2Grid(2,2)

        mapParticle2Grid(1,1) = rp(1)*rp(2)*U*weight*Z
        mapParticle2Grid(1,2) = rp(1)*(1-rp(2))*U*weight*Z
        mapParticle2Grid(2,1) = (1-rp(1))*rp(2)*U*weight*Z
        mapParticle2Grid(2,2) = (1-rp(1))*(1-rp(2))*U*weight*Z
    end function mapParticle2Grid

    function mapGrid2Particle(rp,U)
        implicit none
        real(kind=8) :: rp(2), U(2,2)
        real(kind=8) :: mapGrid2Particle

        mapGrid2Particle = rp(1)*rp(2)*U(1,1)+ &
                           rp(1)*(1-rp(2))*U(1,2)+ &
                           (1-rp(1))*rp(2)*U(2,1)+ &
                           (1-rp(1))*(1-rp(2))*U(2,2)
    end function mapGrid2Particle

    function mapBfield2Grid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), mapBfield2Grid(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_bx(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), temp_by(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: temp_bz(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)

        if(loc(U)==loc(Bz) .or. loc(U)==loc(Bfz) .or. loc(U)==loc(Bwz))then
            ! temp_bz(-ncx:0,-ncy:0) = U(-1-ncx:-1,-1-ncy:-1)
            ! temp_bz(1:Nx+ncx+1,-ncy:0) = U(1:Nx+ncx+1,-1-ncy:-1)
            ! temp_bz(-ncx:0,1:Ny+ncy+1) = U(-1-ncx:-1,1:Ny+ncy+1)
            ! temp_bz(1:Nx+ncx+1,1:Ny+ncy+1) = U(1:Nx+ncx+1,1:Ny+ncy+1)
            temp_bz = getBzGrid(U)
            mapBfield2Grid = (temp_bz(-ncx:Nx+ncx,-ncy:Ny+ncy)+temp_bz(-ncx+1:Nx+ncx+1,-ncy+1:Ny+ncy+1) &
                            + temp_bz(-ncx+1:Nx+ncx+1,-ncy:Ny+ncy)+temp_bz(-ncx:Nx+ncx,-ncy+1:Ny+ncy+1))/4.0d0
        elseif(loc(U)==loc(By) .or. loc(U)==loc(Bfy) .or. loc(U)==loc(Bwy))then
            ! temp_by(-ncx:0,:) = U(-1-ncx:-1,:)
            ! temp_by(1:Nx+ncx+1,:) = U(1:Nx+ncx+1,:)
            temp_by = getByGrid(U)
            mapBfield2Grid = (temp_by(-ncx:Nx+ncx,-ncy:Ny+ncy)+temp_by(-ncx+1:Nx+ncx+1,-ncy:Ny+ncy))/2.0d0
        elseif(loc(U)==loc(Bx) .or. loc(U)==loc(Bfx) .or. loc(U)==loc(Bwx))then
            ! temp_bx(:,-ncy:0) = U(:,-1-ncy:-1)
            ! temp_bx(:,1:Ny+ncy+1) = U(:,1:Ny+ncy+1)
            temp_bx = getBxGrid(U)
            mapBfield2Grid = (temp_bx(-ncx:Nx+ncx,-ncy:Ny+ncy)+temp_bx(-ncx:Nx+ncx,-ncy+1:Ny+ncy+1))/2.0d0
        end if

    end function mapBfield2Grid

    function mapEfield2Grid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), mapEfield2Grid(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_bx(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), temp_by(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
    
        if(loc(U)==loc(Ez))then
            ! write(*,*) "Ez have been on the Grid"
            mapEfield2Grid = U(-ncx:Nx+ncx,-ncy:Ny+ncy)
        elseif(loc(U)==loc(Ex))then
            ! temp_by(-ncx:0,:) = U(-1-ncx:-1,:)
            ! temp_by(1:Nx+ncx+1,:) = U(1:Nx+ncx+1,:)
            temp_by = getByGrid(U)  ! Ex is on the By grid
            mapEfield2Grid = (temp_by(-ncx:Nx+ncx,-ncy:Ny+ncy)+temp_by(-ncx+1:Nx+ncx+1,-ncy:Ny+ncy))/2.0d0
        elseif(loc(U)==loc(Ey))then
            ! temp_bx(:,-ncy:0) = U(:,-1-ncy:-1)
            ! temp_bx(:,1:Ny+ncy+1) = U(:,1:Ny+ncy+1)
            temp_bx = getBxGrid(U)  ! Ey is on the Bx grid
            mapEfield2Grid = (temp_bx(-ncx:Nx+ncx,-ncy:Ny+ncy)+temp_bx(-ncx:Nx+ncx,-ncy+1:Ny+ncy+1))/2.0d0
        end if
    
    end function mapEfield2Grid

    function mapGrid2BxGrid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), mapGrid2BxGrid(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)

        mapGrid2BxGrid = (U(:,-1-ncy:Ny+ncy)+U(:,-ncy:Ny+ncy+1))/2.0d0
    end function mapGrid2BxGrid

    function mapGrid2ByGrid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), mapGrid2ByGrid(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)

        mapGrid2ByGrid = (U(-1-ncx:Nx+ncx,:)+U(-ncx:Nx+ncx+1,:))/2.0d0
    end function mapGrid2ByGrid

    function mapGrid2BzGrid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), mapGrid2BzGrid(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)

        mapGrid2BzGrid = (U(-1-ncx:Nx+ncx,-1-ncy:Ny+ncy)+U(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1) &
        + U(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy)+U(-1-ncx:Nx+ncx,-ncy:Ny+ncy+1))/4.0d0
    end function mapGrid2BzGrid

    function getBxGrid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), getBxGrid(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)

        getBxGrid(:,-ncy:0) = U(:,-1-ncy:-1)
        getBxGrid(:,1:Ny+ncy+1) = U(:,1:Ny+ncy+1)
    end function getBxGrid

    function getByGrid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), getByGrid(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)

        getByGrid(-ncx:0,:) = U(-1-ncx:-1,:)
        getByGrid(1:Nx+ncx+1,:) = U(1:Nx+ncx+1,:)
    end function getByGrid

    function getBzGrid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), getBzGrid(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)

        getBzGrid(-ncx:0,-ncy:0) = U(-1-ncx:-1,-1-ncy:-1)
        getBzGrid(1:Nx+ncx+1,-ncy:0) = U(1:Nx+ncx+1,-1-ncy:-1)
        getBzGrid(-ncx:0,1:Ny+ncy+1) = U(-1-ncx:-1,1:Ny+ncy+1)
        getBzGrid(1:Nx+ncx+1,1:Ny+ncy+1) = U(1:Nx+ncx+1,1:Ny+ncy+1)
    end function getBzGrid

!!!!!!! Generating Ions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    function getMaxwell(A,Ti)
        implicit none
        real(kind=8) :: getMaxwell(3)
        real(kind=8) :: Ti
        integer(kind=4) :: A
        real(kind=8) :: Random(3), Temp

        Ti = Ti/A/938.2d6
        call random_seed()
        call random_number(Random)

        Temp = sqrt(-2.0d0*Ti*log(Random(1)))
        getMaxwell(1) = Temp*cos(2*pi*Random(2))*sin(pi*Random(3))
        getMaxwell(2) = Temp*sin(2*pi*Random(2))*sin(pi*Random(3))
        getMaxwell(3) = Temp*cos(pi*Random(3))
    end function getMaxwell

    function getIonsPosition(i,j)
        implicit none
        real(kind=8) :: getIonsPosition(2)
        integer(kind=4) :: i, j
        real(kind=8) :: Random(2)
        call random_seed()
        call random_number(Random)

        getIonsPosition(1) = dx*Random(1)+x(i,j)
        getIonsPosition(2) = dy*Random(2)+y(i,j)
    end function getIonsPosition

    function cross_product(a,b)
        implicit none
        real(kind=8) :: a(3), b(3)
        real(kind=8) :: cross_product(3)

        cross_product(1) = a(2)*b(3)-a(3)*b(2)
        cross_product(2) = a(3)*b(1)-a(1)*b(3)
        cross_product(3) = a(1)*b(2)-a(2)*b(1)
    end function cross_product

    function cross_product_arr2(ax,ay,az,bx,by,bz,nx,ny)
        implicit none
        real(kind=8) :: ax(nx,ny), bx(nx,ny), ay(nx,ny), by(nx,ny), az(nx,ny), bz(nx,ny)
        real(kind=8) :: cross_product_arr2(nx,ny,3)
        real(kind=8) :: A(3), B(3), C(3)
        integer(kind=4) :: i, j, nx, ny

        do j=1,ny
            do i=1,nx
                A(1) = ax(i,j)
                A(2) = ay(i,j)
                A(3) = az(i,j)
                B(1) = bx(i,j)
                B(2) = by(i,j)
                B(3) = bz(i,j)
                C = cross_product(A,B)
                cross_product_arr2(i,j,1) = C(1)
                cross_product_arr2(i,j,2) = C(2)
                cross_product_arr2(i,j,3) = C(3)
            end do
        end do

    end function cross_product_arr2

    function MAT3_INV(M)
        implicit none
        real(kind=8) :: M(3,3), MAT3_INV(3,3)
        
        MAT3_INV = reshape([M(2,2)*M(3,3)-M(2,3)*M(3,2),M(1,3)*M(3,2)-M(1,2)*M(3,3),M(1,2)*M(2,3)-M(1,3)*M(2,2) &
                           ,M(2,3)*M(3,1)-M(2,1)*M(3,3),M(1,1)*M(3,3)-M(1,3)*M(3,1),M(2,1)*M(1,3)-M(1,3)*M(2,1) &
                           ,M(2,1)*M(3,2)-M(2,2)*M(3,1),M(1,2)*M(3,1)-M(1,1)*M(3,2),M(1,1)*M(2,2)-M(1,2)*M(2,1)],[3,3]) &
        /( M(1,1)*(M(2,2)*M(3,3)-M(2,3)*M(3,2)) - M(2,1)*(M(1,2)*M(3,3)-M(1,3)*M(3,2)) + M(3,1)*(M(1,2)*M(2,3)-M(1,3)*M(2,2)) )

    end function MAT3_INV

    function besselI0(x)    ! polar
        implicit none
        real(kind=8) :: x, y, z, besselI0
        ! if(x==0.0d0) write(*,*) "WARNING! Inf at ZERO"
        if (x<3.75) then
            y = (x/3.75)**2
            besselI0 = 1.0+3.5156229*y+3.0899424*y**2+1.2067492*y**3 &
            +0.2659732*y**4+0.0360768*y**5+0.0045813*y**6
        else
            z = 3.75/x
            besselI0 = exp(x)*(0.39894228+0.01328592*z+0.00225319*z**2 &
            -0.00157565*z**3+0.00916281*z**4-0.02057706*z**5 &
            +0.0263557*z**6-0.01647633*z**7+0.00392377*z**8)/sqrt(x)
        end if


    end function besselI0

    function besselk0(x)    ! polar
        implicit none
        real(kind=8) :: x, y, z, besselk0
        ! if(x==0.0d0) write(*,*) "WARNING! Inf at ZERO"
        if (x<2) then
            y = (x/2.0)**2
            besselk0 = -0.57721566 + 0.4227842*y + 0.23069756*y**2 + 0.0348859*y**3 &
            + 0.00262698*y**4 + 0.0001075*y**5 + 0.0000074*y**6-besselI0(x)*log(0.5*x)
        else
            z = 2.0/x
            besselk0 = exp(-x)*(1.25331414-0.07832358*z+0.02189568*z**2 &
            -0.01062446*z**3+0.00587872*z**4-0.0025154*z**5 +0.00053208*z**6)/sqrt(x)
        end if

    end function besselk0

!!!!!!! Recording Elapsed Time !!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine recordStartTime()
        implicit none

        StartTime = mpi_wtime()
    end subroutine recordStartTime

    subroutine recordEndTime()
        implicit none

        call mpi_barrier(mpi_comm_world,ierr)
        EndTime = mpi_wtime()
    end subroutine recordEndTime

    subroutine displayUsedTime()
        implicit none

        if (myid==0) then
            write(*,*) "Elapsed time is ", EndTime-StartTime ," seconds."
        end if
    end subroutine displayUsedTime

!!!!!!! Release Memory !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine releaseMemory()
        implicit none
        type(PseudoParticle), pointer :: p

        do while(associated(IonListHead))
            p => IonListHead
            IonListHead => IonListHead%next
            deallocate(p)
        end do
        
        do while(associated(InjectorListHead))
            p => InjectorListHead
            InjectorListHead => InjectorListHead%next
            deallocate(p)
        end do

        do while(associated(MPIListHeadLeft))
            p => MPIListHeadLeft
            MPIListHeadLeft => MPIListHeadLeft%next
            deallocate(p)
        end do

        do while(associated(MPIListHeadRight))
            p => MPIListHeadRight
            MPIListHeadRight => MPIListHeadRight%next
            deallocate(p)
        end do

        do while(associated(MPIListHeadUp))
            p => MPIListHeadUp
            MPIListHeadUp => MPIListHeadUp%next
            deallocate(p)
        end do

        do while(associated(MPIListHeadDown))
            p => MPIListHeadDown
            MPIListHeadDown => MPIListHeadDown%next
            deallocate(p)
        end do

        do while(associated(MPIInjectorListHeadLeft))
            p => MPIInjectorListHeadLeft
            MPIInjectorListHeadLeft => MPIInjectorListHeadLeft%next
            deallocate(p)
        end do

        do while(associated(MPIInjectorListHeadRight))
            p => MPIInjectorListHeadRight
            MPIInjectorListHeadRight => MPIInjectorListHeadRight%next
            deallocate(p)
        end do

        do while(associated(MPIInjectorListHeadUp))
            p => MPIInjectorListHeadUp
            MPIInjectorListHeadUp => MPIInjectorListHeadUp%next
            deallocate(p)
        end do

        do while(associated(MPIInjectorListHeadDown))
            p => MPIInjectorListHeadDown
            MPIInjectorListHeadDown => MPIInjectorListHeadDown%next
            deallocate(p)
        end do

        ! Deallocate Grid Variables
        if(allocated(x)) deallocate(x)
        if(allocated(y)) deallocate(y)
        if(allocated(Ex)) deallocate(Ex)
        if(allocated(Ey)) deallocate(Ey)
        if(allocated(Ez)) deallocate(Ez)
        if(allocated(Exg)) deallocate(Exg)
        if(allocated(Eyg)) deallocate(Eyg)
        if(allocated(Bx0)) deallocate(Bx0)
        if(allocated(By0)) deallocate(By0)
        if(allocated(Bz0)) deallocate(Bz0)
        if(allocated(Bx)) deallocate(Bx)
        if(allocated(By)) deallocate(By)
        if(allocated(Bz)) deallocate(Bz)
        if(allocated(Bxg)) deallocate(Bxg)
        if(allocated(Byg)) deallocate(Byg)
        if(allocated(Bzg)) deallocate(Bzg)
        if(allocated(Bwx)) deallocate(Bwx)
        if(allocated(Bwy)) deallocate(Bwy)
        if(allocated(Bwz)) deallocate(Bwz)
        if(allocated(Bfx)) deallocate(Bfx)
        if(allocated(Bfy)) deallocate(Bfy)
        if(allocated(Bfz)) deallocate(Bfz)
        if(allocated(ux)) deallocate(ux)
        if(allocated(uy)) deallocate(uy)
        if(allocated(uz)) deallocate(uz)
        if(allocated(ne)) deallocate(ne)
        if(allocated(ni)) deallocate(ni)
        if(allocated(ne_prev)) deallocate(ne_prev)
        if(allocated(ne_temp)) deallocate(ne_temp)
        if(allocated(Te)) deallocate(Te)
        if(allocated(pe)) deallocate(pe)
        if(allocated(ne0)) deallocate(ne0)
        if(allocated(Te0)) deallocate(Te0)
        if(allocated(pe0)) deallocate(pe0)
    end subroutine releaseMemory
end module Setup
